#include"mrrr.h"
int pdlaneg(int *n, double *d__, double *lld, double * sigma, double *pivmin,
		int *r__) {
	/* System generated locals */
	int ret_val, i__1, i__2, i__3, i__4;

	/* Local variables */
	int j;
	double p, t;
	int bj;
	double tmp;
	int neg1, neg2;
	double bsav, gamma, dplus;
	int negcnt;
	bool sawnan;
	double dminus;

	/*  -- LAPACK auxiliary routine (version 3.2) -- */
	/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
	/*     November 2006 */

	/*     .. Scalar Arguments .. */
	/*     .. */
	/*     .. Array Arguments .. */
	/*     .. */

	/*  Purpose */
	/*  ======= */

	/*  DLANEG computes the Sturm count, the number of negative pivots */
	/*  encountered while factoring tridiagonal T - sigma I = L D L^T. */
	/*  This implementation works directly on the factors without forming */
	/*  the tridiagonal matrix T.  The Sturm count is also the number of */
	/*  eigenvalues of T less than sigma. */

	/*  This routine is called from DLARRB. */

	/*  The current routine does not use the PIVMIN parameter but rather */
	/*  requires IEEE-754 propagation of Infinities and NaNs.  This */
	/*  routine also has no input range restrictions but does require */
	/*  default exception handling such that x/0 produces Inf when x is */
	/*  non-zero, and Inf/Inf produces NaN.  For more information, see: */

	/*    Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in */
	/*    Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on */
	/*    Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624 */
	/*    (Tech report version in LAWN 172 with the same title.) */

	/*  Arguments */
	/*  ========= */

	/*  N       (input) int */
	/*          The order of the matrix. */

	/*  D       (input) DOUBLE PRECISION array, dimension (N) */
	/*          The N diagonal elements of the diagonal matrix D. */

	/*  LLD     (input) DOUBLE PRECISION array, dimension (N-1) */
	/*          The (N-1) elements L(i)*L(i)*D(i). */

	/*  SIGMA   (input) DOUBLE PRECISION */
	/*          Shift amount in T - sigma I = L D L^T. */

	/*  PIVMIN  (input) DOUBLE PRECISION */
	/*          The minimum pivot in the Sturm sequence.  May be used */
	/*          when zero pivots are encountered on non-IEEE-754 */
	/*          architectures. */

	/*  R       (input) int */
	/*          The twist index for the twisted factorization that is used */
	/*          for the negcount. */

	/*  Further Details */
	/*  =============== */

	/*  Based on contributions by */
	/*     Osni Marques, LBNL/NERSC, USA */
	/*     Christof Voemel, University of California, Berkeley, USA */
	/*     Jason Riedy, University of California, Berkeley, USA */

	/*  ===================================================================== */

	/*     .. Parameters .. */
	/*     Some architectures propagate Infinities and NaNs very slowly, so */
	/*     the code computes counts in BLKLEN chunks.  Then a NaN can */
	/*     propagate at most BLKLEN columns before being detected.  This is */
	/*     not a general tuning parameter; it needs only to be just large */
	/*     enough that the overhead is tiny in common cases. */
	/*     .. */
	/*     .. Local Scalars .. */
	/*     .. */
	/*     .. Intrinsic Functions .. */
	/*     .. */
	/*     .. External Functions .. */
	/*     .. */
	/*     .. Executable Statements .. */
	/* Parameter adjustments */
	--lld;
	--d__;

	/* Function Body */
	negcnt = 0;
	/*     I) upper part: L D L^T - SIGMA I = L+ D+ L+^T */
	t = -(*sigma);
	i__1 = *r__ - 1;
	for (bj = 1; bj <= i__1; bj += 128) {
		neg1 = 0;
		bsav = t;
		/* Computing MIN */
		i__3 = bj + 127, i__4 = *r__ - 1;
		i__2 = min(i__3, i__4);
		for (j = bj; j <= i__2; ++j) {
			dplus = d__[j] + t;
			if (dplus < 0.) {
				++neg1;
			}
			tmp = t / dplus;
			t = tmp * lld[j] - *sigma;
			/* L21: */
		}
		sawnan = isnan(t);
		/*     Run a slower version of the above loop if a NaN is detected. */
		/*     A NaN should occur only with a zero pivot after an infinite */
		/*     pivot.  In that case, substituting 1 for T/DPLUS is the */
		/*     correct limit. */
		if (sawnan) {
			neg1 = 0;
			t = bsav;
			/* Computing MIN */
			i__3 = bj + 127, i__4 = *r__ - 1;
			i__2 = min(i__3, i__4);
			for (j = bj; j <= i__2; ++j) {
				dplus = d__[j] + t;
				if (dplus < 0.) {
					++neg1;
				}
				tmp = t / dplus;
				if (isnan(tmp)) {
					tmp = 1.;
				}
				t = tmp * lld[j] - *sigma;
				/* L22: */
			}
		}
		negcnt += neg1;
		/* L210: */
	}

	/*     II) lower part: L D L^T - SIGMA I = U- D- U-^T */
	p = d__[*n] - *sigma;
	i__1 = *r__;
	for (bj = *n - 1; bj >= i__1; bj += -128) {
		neg2 = 0;
		bsav = p;
		/* Computing MAX */
		i__3 = bj - 127;
		i__2 = max(i__3, *r__);
		for (j = bj; j >= i__2; --j) {
			dminus = lld[j] + p;
			if (dminus < 0.) {
				++neg2;
			}
			tmp = p / dminus;
			p = tmp * d__[j] - *sigma;
			/* L23: */
		}
		sawnan = isnan(p);
		/*     As above, run a slower version that substitutes 1 for Inf/Inf. */

		if (sawnan) {
			neg2 = 0;
			p = bsav;
			/* Computing MAX */
			i__3 = bj - 127;
			i__2 = max(i__3, *r__);
			for (j = bj; j >= i__2; --j) {
				dminus = lld[j] + p;
				if (dminus < 0.) {
					++neg2;
				}
				tmp = p / dminus;
				if (isnan(tmp)) {
					tmp = 1.;
				}
				p = tmp * d__[j] - *sigma;
				/* L24: */
			}
		}
		negcnt += neg2;
		/* L230: */
	}

	/*     III) Twist index */
	/*       T was shifted by SIGMA initially. */
	gamma = t + *sigma + p;
	if (gamma < 0.) {
		++negcnt;
	}
	ret_val = negcnt;
	return ret_val;
} /* dlaneg_ */
